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ABSTRACT 

We present a scaling law that predicts the values of the stresses obtained in numerical simulations of saturated 
MRI-driven turbulence in non-stratified shearing boxes. It relates the turbulent stresses to the strength of the 
vertical magnetic field, the sound speed, the vertical size of the box, and the numerical resolution and predicts 
accurately the results of 35 numerical simulations performed for a wide variety of physical conditions. We use 
our result to show that the saturated stresses in simulations with zero net magnetic flux depend linearly on the 
numerical resolution and would become negligible if the resolution were set equal to the natural dissipation 
scale in astrophysical disks. We conclude that, in order for MRI-driven turbulent angular momentum transport 
to be able to account for the large value of the effective alpha viscosity inferred observationally, the disk 
must be threaded by a significant vertical magnetic field and the turbulent magnetic energy must be in near 
equipartition with the thermal energy. This result has important implications for the spectra of accretion disks 
and their stability. 

Subject headings: black hole physics — accretion, accretion disks — MHD — instability — turbulence 



1. INTRODUCTION 

One of the key unsolved problems in accretion disk physics 
is the precise nature of the mechanism that allows for an- 
gular momentum to be transported outwards in the disk. In 
the standard theory, ones copes with this problem by argu- 
ing that the stress responsible for angular momentum trans- 
port, due to turbulence and magnetic fields, is proportional 
to the local pressure, i.e., T r< p = qaP. Here, a is a 
dimensionless constant of the order of, but smaller than, 
unity and q = — dlnfl/dlnr characterizes the local shear 
( Shak ura & Sunvaevlll973HKato. Fukue. & Mineshigdll998l: 
iFrank. King. & Rainell2002h . 

Despite the fact that the standard parametrization leads to a 
disk model in which the energy generation rate is determined 
mostly by en ergy balance and depends we akly on the adopted 
prescription (Balbus & Papal oizoul 1 19991) . almost every as- 
pect of the disk structure depends explicitly on this assump- 
tion (j Kato. Fukue, & Mineshige 1998; Frank, King, & Raind 
2002). Therefore calculating the value of a, or more gener- 
ally, testing whether the assumed relationship between stress 
and pressure is adequate, is of fundamental importance. These 
questions, however, lie outside the scope of the standard the- 
ory, making necessary the identification of a specific physical 
mechanism for angular momentum transport. 

Over the last decade, magnetohydrodynamic (MHD) tur- 
bulence driven by the magnetorotational instability (MRI; 
Balbus & Hawley 1991; 1998) has emerged as the most 
promising candidate to enable angular momentum trans- 
port in astrophysical disks. The development of three- 
dimensional MHD numerical codes has led to a detailed 
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study and characterization of angular momentum transport 
in turbulent magnetized disks. This has motivated nu- 
merical estimations of the a parameter and, more gen- 
erally, the search for saturatio n predictors to d e scribe 
the turbulent state (see, e.g., lHawlev etafl 1 1 9951 119961; 
i Brandenburg. Nordlund. Stein. & Torkelssonlll995l ; iGammie 
1998;lPessah. Chan. & Psalfisll2007l) . " 

It has long been recognized that the a parameter is not 
a constant. It is known to depend, among other things, on 
the strength and geometry of the magnetic field and, perhaps 
more uncomfortably, on the size of the simulation domain and 
the resolution. In spite of this, it has not been possible to 
disentangle nu merical from physical dependencies in a clear 
way (see, e.g.. lHawlev et alii 19951 1 1996b iBrandenburdl 1998b 
ISano. Inutsuka. & Mivamal Il998t ISano et all 120041) . Being 
able to distinguish between these two types of dependen- 
cies is vital if we seek to use the results of numerical sim- 
ulations to build angular momentum transport models, and 
eventually global d isk models, beyond the standard pre- 
scription (see, e.g.. iKato & Yoshizawalll995l : lOgilvid [20031 : 
iPessah. Chan. & Psaltisll2006bll2007l) . 

The search for a mechanism for the saturation of the MRI- 
driven turbulence has been a long-sought for goal since the ap- 
preciation of the relevance of the MRI to astrophysical disks. 
By necessity, a key piece of this puzzle consists of understand- 
ing the role that the various factors, both physical and numer- 
ical, play in the saturated state. In this Letter, we provide an 
expression that describes the transport of angular momentum 
in shearing MHD boxes, b ased on a se r ies of local numerical 
simulations carried out by ISano et al.l (120041) . In particular, 
we are able to disentangle how the characteristics of the MHD 
flow depend on the various physical (pressure, magnetic field, 
etc.) and numerical (box size and resolution) factors. 

2. CHARACTERISTIC SCALES IN NUMERICAL SIMULATIONS 

Numerical studies addressing the local dynamics of three- 
dimensional, differentially rotating, turbulent magnetized 
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flows are often carried out in the shearing box approximation. 
This consists of a first order expansion in the variable r — ro 
of all the quantities characterizing the flow at the fiducial ra- 
dius tq. The goal of this approach is to retain the most relevant 
physics governing the dynamics of the MHD fluid in a locally- 
Cartesian coordinate system co-orbiting and corotating with 
the background flow with local velocity Vq = tq f2o0. 

In the shearing box framework, all the physical variables 
are usually normalized using the initial density, po, and the 
angular velocity, fio- For a more detailed discussion concern- 
ing the physical approximations and numerica l implementa- 
tions i nvolved in the shearing box approach see lHawley et al.1 
( 1995). In an unstratified shearing box there are three possi- 
ble scales of length relevant to the vertical direction: ( i) The 
first scale is the size of the box L. ( ii) The second scale is the 
wavelength corresponding to the most unstable MRI mode, 

i.e., Amri = 27r-y/l6/15uAz/^0i where vaz is the average 
Alfven velocity associated with the vertical component of the 
magnetic field, vaz = B z /(4 : ttpq) 1 ^ 2 , and we have assumed 
a Keplerian shear profile (i.e., q — 3/2). (Hi) The last scale is 
associated with the sound speed: H = (2/j) 1 ' 2 c s /Qo, where 
7 is the ratio of specific heats and c s is the speed of sound. In 
a stratified disk, H is the pressure scale-height and, for sim- 
plicity, we will refer to it as such hereafter, even though the 
particular simulations we will be discussing are not stratified. 

Most systematic studies carried out to characterize the tur- 
bulent state driven by the MRI consider non-radiative flows in 
which the internal energy grows in time due to magnetic dissi- 
pation. Furthermore, most of the initial efforts to characterize 
the saturated state of the MRI employed numerical schemes 
that evolved the internal energy, as opposed to the total en- 
ergy. In these cases the reported values of a, as well as the 
expressions for different predictor functions, where obtained 
considering the initial pressure Po- In such simulations, how- 
ever, gas pressure increases linearly with time and this effect 
must be considered in order to understand its effects on the 
stresses responsible for angular momentum transport. 

ISano et aT] ((2004), carried out a series of numerical simu- 
lations in order to investigate the effects of the evolving pres- 
sure on the turbulent state. They employed an algorithm that 
solves the energy equation in terms of the total energy, which 
allowed them to keep better track of the energy budget of the 
flow. By considering the vertical extent of the box L as the 
unit of length, they were able to examine the dependence of 
the efficiency of angular momentum transport on the gas pres- 
sure. In this case, both the gas pressure Po and the mean 
Alfven velocity va z , are independent parameters and deter- 
mine the ratios Hq/L and Amri/P respectively, where the 
subscript zero indicates initial values. 

3. SCALING LAWS IN MRI-DRIVEN TURBULENCE 



ISano et ail d2004l) performed an extensive set of simula- 
tions with zero net magnetic flux through the domain, initial- 
ized with a vertical magnetic field Bo z (r) = Bq z sin[27r(r — 
ro)/L r ], and runs with a uniform magnetic field, Bq z , per- 
pendicular to the disk midplane 5 . They investigated adiabatic 
(7 = 5/3) as well as isothermal (7 = 1.001) equations of 
state. They also performed a suite of runs to study the ef- 
fects of Ohmic dissipation but we will not discuss them here. 

5 Note that in both of these cases the magnetic flux through the domain is 
conserved because of the adopted vertical periodic boundary conditions. 
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FIG. 1. — Dimensionless magnetic stress, normalized by the instantaneous 
pressure, as a function of a dimensionless measure of the box size with re- 
spect to the ev olving scale-heigh t, H = Hq(P / Pq) 1 / 2 , for the numerical 
simulations of Sano et al. (2004). The various types of symbols label sim- 
ulations according to the value of the ratio v\ z / LQo, with L = 1 and 
On = 10 -3 . The dotted lines are the best fits for each set of simulations 
with constant ratio v\ z / LSlo 7^ 0. The dashed line is the best fit for all the 
simulations with zero net magnetic flux. In all the cases a constant slope of 
5/3, i.e., Mj.^/ P oc (L/H) 5 / 3 , provides a remarkably good description of 
the simulation results. 

For all their models they used a shearing box with dimen- 
sions L = L r — L z = 1 and L$ = 4X, and a grid of 
32 x 128 x 32 zones. The sca les of density and time are the 
same as in lHawlev et al.l (11995b . i.e., p Q — 1 and fl a = 10~ 3 . 

Their choice of initial conditions spans six orders of mag- 
nitude in the initial pressure Po and two orders of magnitude 
in the Alfven speed. Figure [T] shows the dimensionless ra- 
tios M rt f,/P as a function of the ratio L/H characterizing the 
turbulent states reached by the adia batic and isotherm al sim- 
ulations listed in Tables l and 2 in ISano et all (120041) . Here, 
M r = ((5B r 5Bj,))/4:ir, P, and H stand for the volume- and 
time- averaged values of the magnetic stress, pressure, and, 
equivalent scale-height H = H (P/Po) 1 / 2 . 

The tight correlations followed by simulations character- 
ized by the same values of mean Alfven velocities, includ- 
ing the class of simulations with zero net magnetic flux, i.e., 
v Az = 0, suggest the scaling 
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The best fits to the data are shown in FigureQ]with a dashed 
line and dotted lines for simulations with zero and non-zero 
net magnetic flux, respectively. This scaling 6 can be used to 
remove the dependence of the ratio M r< f, / P on L/H for each 
set of simulations with the same value of the ratio vaz/L^o- 
The result is shown in Figure [2] where the various types of 
symbols label the best fit values characterizing each class of 
simulations according to their value of vazI 'LQ,q, The error 
bars quantify the scatter within each class of simulations from 
the best fit values obtained in Figure Q] These average values 
follow a simple, yet tight, correlation with the associated val- 
ues of Amri- 

6 Note that fixing all the slopes to 5/3 in Fig.[T]allows us to describe all the 
simulations on the same footing while leading to small residuals in the best 
fit amplitudes (see Fig. [5} 
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FlG. 2. — Dimensionless magnetic stress, —M r< f,/P, multiplied by 
(H I X) 5 / 3 , as a function of the wavelength corresponding to the most un- 
stable MRI mode, Amri- The various types of symbols correspond to the 
best fit values characterizing each class of simulations according to the value 
of the ratio v\ z /LQo, as inferred from FigureffJ The simulations with zero 
mean magnetic flux, i.e., vaz = 0, are displayed at some arbitrary value for 
visualization purposes only. The error bars quantify the scatter within each 
class of simulations around the corresponding mean values. Vertical dotted 
lines represent the values at which the most unstable MRI wavelength equals 
the grid size and the size of the box, respectively, i.e., Amri = A = 1/32 
and Amri = L = 1. The overall dependence of the saturation level on the 
ratio Amri/£ (solid line) is given by the saturation predictor ^3). Numerical 
simulations for which Amri > L are stable to the MRI (shaded region). 

The vertical dotted lines in Figure [2]represent the values at 
which the most unstable MRI wavelength equals the size of 
the box and the grid size, respectively, i.e., Amri = L = 1 
and Amri = A = 1/32. Note that numerical simulations for 
which Amri > L are stable to the MRI (shaded region). The 



saturation of the simulations for which A < 
linearly proportional to Amri, i-e-, 
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This steeper dependence does indeed provide a good 



description of the simulations with va z = 5 x 10~ 5 ,2.5 x 
10~ 5 , 1.25 x 10~ 5 . However, a linear dependence on va z 
leads to a better overall description of all the simulations with 
A < Amri < L, i.e., including those with vaz = 6.25 x 10~ 6 
and vaz = 1 x 10~ 4 . 

There is a clear departure from the linear dependence of 
the saturated stress on va z for the simulations for which 
A 



MRI 



< A. All of these runs saturate at the same value of 
{M r(t ,/P)(H/L f/' i regardless of the value of the magnetic 
flux through the vertical boundary. The fact that this mini- 
mum value of the stresses at saturation (lower dashed line in 
Figure [2]i is equal to 1/32 of the maximum possible value, 
corresponding to Amri = L (upper dashed line in the same 
Figure), strongly suggests that this floor is entirely set by the 
grid size. Extrapolating this behavior to lower values of the 
grid scale suggests that this minimum saturation level is itself 
linearly proportional to the size of the grid 7 . This is consistent 

7 After this paper appeared on the preprint server, Fromang and Papaloizou 
posted a paper I arXiv:0705.3621 1 in which they perform ideal MHD simula- 
tions for zero-net-field shearing boxes with increasing resolution. They found 



with the results found in lHawlev et alj (|T996); see in particu- 
lar their Figure 8, where the final plasma /3's for simulations 
with zero net magnetic flux and uniform vertical magnetic 
fields is shown. The ratio of the final plasma (3 for the sim- 
ulation with highest uniform field to the average final plasma 
j3 corresponding to the zero net magnetic flux runs is indeed 
close to A = 1/32. 

The overall dependence of the saturation level on the ratio 
Amri/^, i-e., the solid line in Figure |2] is described by the 
function 
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This saturation predictor is consi stent, in the interm ediate 
regime, with the one obtained by Hawley et al.l (Il995h that 
we also used in an earlier work (Pessah et al. 2006b). In- 
deed, substituting in this expression P = Pq = po{HqVL^) 2 , 
H = Hq = L, and considering that the magnetic stress 
is_ roughly half of t he magnetic energy in the fluctu ations, 
M = {{5B 2 ))/8tt dBlackman. Penna. & Varnierell2006l and 
references therein), leads to 

M = 1.2p (i«o) (Amri^o) , (4) 

which is identical to equation (18) of lHawley et al.l (1 19951) . 

The 5/3 scaling in the saturation predictor (fJJ is reminis- 
cent of the Kolmogorov spectrum of turbulence. The promi- 
nent role of Amri in this expression is another indication that 
the MRI continues to pump the turbulence, even in the non- 
linear state (see also Pessah et al. 2006a, 2006b). The depen- 
dence of the stress on the vertical size of the box, L, must 
be an artifact of the periodic boundary conditions, which do 
not allow for any turbulent energy to escape the domain of 
solution. The dependence on the resolution must be related 
to numerical dissipation at the grid scale, and might, there- 
fore, change if a different MHD algorithm were used for the 
simulations. Finally, the role of the "scale-height", H, is very 
difficult to understand, since the simulations are not stratified 
and this length scale is too large to be resolved. It may arise 
from a non-linear coupling between sound waves and MHD 
modes or it might simply be related to numerical dissipation 
through the dependence of the Courant time step on sound 
speed. 

4. IMPLICATIONS AND DISCUSSION 

In shearing-box simulations of Keplerian flows, the 
Maxwell and Reynol ds, R r( t, = ((pSvrSvd,)), stresses follow a 
tight correlation (see Pessah, Chan. & P saltis 2006a, and ref- 
erences therein) with —M r ^l R r< §, ~ 4. Taking this into ac- 
count, we can write for the total stress and the effective alpha 
viscosity 

A/L if Amri < A 
Amri/ L if A < Amri < L , 
if A M ri > L 

(5) 

where q = 3/2 for a Keplerian disk. This is a remarkable 
result. This expression accurately describes the overall de- 
pendence of the saturated state for 35 numerical simulations 

that the saturated stresses depend linearly on resolution, in agreement with the 
results presented here. 
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spanning six orders of magnitude in initial pressure, encom- 
passing domains with zero and non-zero net magnetic flux, as 
well as adiabatic and isothermal equations of state. 

Note that the very small values of the effective alpha vis- 
cosity reported in the past, which ar e unable to account for 
obser vations of astrophysical disks dKing. Pringle. & Livid 
120071) . correspond to shearing box simulations with vertical 
extents that are small compared to the equivalent pressure 
scale-height. One would expect that a realistic value of the 
stress is achieved in simulations with L ~ H. Nevertheless, 
even having taken this into account, large values of the alpha 
viscosity can be achieved only for particular configurations. 
This has three important implications for accretion disk mod- 
els in which the angular momentum transport is mediated by 
MRI-driven turbulence. 

First, the saturated stresses in simulations with zero net 
magnetic flux are linearly proportional to the numerical res- 
olution. This implies that if we were able to set the numeri- 
cal resolution to the natural dissipation scale in the problem, 
which is many orders of magnitude smaller than the pressure 
scale height, the MRI would be unable to sustain the nec- 
essary turbulent stresses in a Keplerian shearing box, unless 
there were a significant magnetic flux through the vertical box 
boundaries. 

Second, in order for MRI-driven turbulence to account for 
the large values of the effective alpha viscosity inferred from 
observations (a > 0.1; see King et al. 2007), the vertical 
magnetic field must grow to a strength such that the most un- 
stable MRI-modes have wavelengths comparable to the disk 
scale-height, i.e., Amri — L ~ H. Such vertical fields 
are only a small fraction, roughly a few hundredths, of the 
associated equipartition field and pose no significant prob- 
lem to the energy budget of the accretion flow. Indeed, the 
origin of such a small mean magnetic field perpendicular to 
the disk midplane need not be external to the disk. MRI- 
driven fluctuations can easily give rise to magnetic fields of 



this order, perhaps through the combined effects of shear- 
ing, MR I, and Parker instabilit y, as in the mechanism pro- 
posed by Tout & Pringle dl992l) . Note that disk stratification 
is likely to play an important role in driving helical turbulence 
and thus in ena bling the development of a global, large-scale 
magnetic flux (|Brandenburg ; Nordlund. Stein. & Torkelssonl 
1995; Tan & Blackman||200ll 

Black man & Tanll2004h . 

Finally, in MRI-driven turb ulence, the turbulent magnetic 
energy is comparable to T rt j, dPessah. Chan. & Psaltis 2006a) 
and, hence, will have to be also comparable to the thermal 
energy, if a is of order unity. This implies that the vertical 
scale height of an accretion disk is set by both the magnetic 
and the thermal pressures, with important i mplication for the 
spectrum emerging from each disk annulus dBlaes et al.l2006l) 
as well as for the viscous and thermal stability of the disk. 

In closing, it is important to remark that the satura- 
tion predictor <(5j may be specific to shearing-box simula- 
tions with periodic boundary conditions and not necessar- 
ily applicable to the local saturation of stresses in strati- 

I II I 

fled (e .g., Brandenburg, Nordlund, Stein . & Torkelssonll 1995c 
Miller & Stone 2000) or global simulations (e.g., Armitage 
1998; Hawlev 2000, 2001). In these cases, the physical mech- 
anism that limits the growth of turbulent magnetic energy 
maybe related to magnetic buoyancy or to large meridional 
circulation. In both mechanisms, magnetic energy is lost at 
large scales and, therefore, the dependence of the saturated 
stress on the numerical resolution and the sound speed may 
disappear. Understanding the physical origin of the satura- 
tion predictor (0 and comparing it to simulations of stratified 
shearing boxes will resolve these issues. 



We thank Eric Blackman, Jeremy Goodman, Gordon 
Ogilvie, and Jim Stone for valuable comments and discus- 
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